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i. eel ERODUC TION 


Several recent attempts to understand the complex 
dynamic interaction of the thunderstorm and its environment 
have been made using the dynamical thunderstorm model pre- 
sented by Newton and Newton (1959). The model proposes that 
Zaeeomage vertical drafts within the thunderstorm create an 
Serective barrier to the environmental flow. The analysis 
by Browning (1964) concerning airflow near and within severe 
thunderstorms provides valuable additional knowledge about 
Ene character of the vertical motion. Convective storm mass 
and water budget assessments such as the one conducted by 
Newton (1966) substantiate the concept of upper level 
“HeeamnlentsOor COOl, dry air first recognized by Normand 
(1946). 

Radar determined chaff trajectories were used success- 
fully by Fankhauser (1968) to follow the motion of air 
parcels in and around a thunderstorm. Further use of chaff 
will undoubtedly make possible a good understanding of the 
paths taken by environmental parcels from all positions 
around a local storm. 

As noted by Lilly (1962), the evolution of convective 
motions with relatively large amplitudes can best be 
described through the combined application of conventional 
analytic methods and numerical experimentation. Ogura (1963) 
discussed the possible use of the primitive equations as a 


basis for forecasting small scale dynamical features. He 


IEE 





suggested that interaction between convective elements and 
the prevailing wind field could thus be investigated numeri- 
Gally. The results of such efforts, correlated with the 
subjective modeling, should produce verifying as well as 
augmenting information for the subjective products. 

A numerical study of turbulent flow conducted by 
Deardorff (1970) makes use of a three dimensional primitive 
€quation model. The forecast procedure he uses is similar 
to the procedure used in this study. 

The numerical model presented here is a three dimen- 
Sional primitive equation model based on the Boussinesq 
e@quations. Forecasts are made for the wind components and 
potential temperature after impoSing initial conditions on 
these same parameters. 

Piewgeiemased £Or Computations is 25 by 25 in the 
horizontal, with 15 levels in the vertical. Grid interval 
in all directions is 1 kilometer. A mean latitude of 35° 
is used to determine the Coriolis parameter, which is 
treated as a constant. 

Predictions are made using the IBM OS/360 computer of 
the W. R. Church computer center. A one minute time step 
was chosen to conform to computational stability require- 


ments. 
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Peay ELOEMENT OF THE MODEL 


A. THE FORECAST EQUATIONS 
the pressure term in the equation of motion may be 


expressed as 





Pp 1/1 
a PRE O 
Pp Ve P= > V3 1/x B 
oP 
=" ©V38 . 
where 
kK 
Pp 
B =C — 
p B. 
and 


@ 
= 
r-— 
rio" 
aa 


Using this result in the component equations of motion 


yields 
2 
du. 9 28 + fv + YV3 4, (1) 
Z 
a = -0 38 ea ee Ce (2) 
wn - 0 28 -g4 yv,? w. (3) 
Z 


Forming the perturbation expressions 


B= B,(2) + B'(x.y,z,t) 


and 


iS 





Ce 0 ot 8 O (x,y,z, t) 


where 
|B" << Bo and |a|<<@, - 


and substituting into (1) yields 


; Z 


’ 
Peece the relatively small product 9 28 is neglected. 
cimatlarly, 


; 2 
dv =~ 9, OF - fu + UV3 Ne (5) 


Mem@ethne vertical component equation the substitution yields 





rae. obo. SL SE 2 
7 uP Aus 0, oon i Sie Te NEG) Ww. 


Since B. and Us are related by the hydrostatic equation 





the vertical component equation may be written 


She op , 28 , ow. 5) 
de 9052 * @ Oe ai 


Defining op = 0, B' and substituting into (4), (5) and 


(6) results in the system of equations 


14 








ou | ap 2 

Se = sett) aS + fv + A ale U (7) 

2) : 

a = iy) dy HEY ap VV, Vv (8) 

g@ 2 
Ow = — (w) - oP +—— + YVV W (9) 
Ont on 3 
Oo 

a 8) (10) 

ait ; 

Vv, : V3 = 0 (11) 
where ¢(S) = “V3° V3 Sees, vw 1G ee cua be lLonm (1) ) sis tne 
SeneiInulty Equation. In order to eliminate sound waves, 


incompressible flow is assumed. Ogura and Phillips (1962) 
have shown that this system of equations, which they call 
Che anelastic equations, is quite accurate for shallow con- 


vection. 


Bs DETERMINATION OF THE PRESSURE TERMS 
Equations (7), (8), and (9) may be combined and written 
in the vector form 


OY 


Taking the three dimensional divergence of (12) yields 


> 








V2°- V 
5 v2 2), Es, 
ot 
meat — De 
GaN tae | (V3 > v3) v3 | - V3 D 
(13) 

; ss may g 06 
Vg | £K . v3) |+ 6, Oz 


For a limited horizontal scale the Coriolis parameter may be 


treated as constant. Thus, (13) becomes 








2 = - “oS G 930 
V; = V3 ° | (V3 * ¥3)V3]- fo V3 ° OK x V5) + Q, 82 
7 _— . — jh a8 
= - V3 : | (V3 V3) V3 4 i Cay 6. Sa) e 1 s) 


Ponca Deusolved using relaxation techniques. 
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Met Lut he DIP FERENCING 


Equations (7) through (10) are solved numerically by 
introducing finite differences in x, y, Z and t. Solution 
for the pressure term can be carried out by applying the 
extrapolated Liebmann relaxation technique to evaluate (14). 
Appendix A shows the development of the relaxation equation. 
The peecedure wsed to compute the forcing function is 
included in Appendix B. 

In order to maintain finite differencing consistency 
throughout the entire forecast process, it is necessary to 
use the finite difference forms of equations (7), (8) and 
(9) to form the finite difference equivalent of equation 
(14). 

The vertical domain consists of fifteen levels with 
Values of u, v, w, @, and @ assigned to each grid point. 
The space differencing method used is based on a scheme 
devised by Arakawa (1966) and is designed to conserve total 
EHergy- Ihe non-linear operator, £(S), takes' the finite 


difference form 
= * 1 + S. ‘ ) 
oe [isa get Mage) Sten ie * Si, 5, 


S 


AS V4 Ax 


; ; : =< we 
Oe eee a ute en 3 eel Dees (or < shag Bre 


miele Val gi) Si jclk Si 5.e)| ao 


ny 





“ Hs jae. t Ma, 5,20) (Oi pp inte 
(ieee * 5.x) (i,5.44 + Si 3,0) /4 Az. 


The linear space derivatives are of the centered form 


ss _ itt, 5% 7 Si-2, i.) 
Ox 2 Ax ; 


Centered time differencing is used for all forecasts except 


Paeerinst which employs a forward time step. 
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IV. BOUNDARY CONDITIONS 


The boundary conditions are selected to enforce mass and 
energy conservation within the domain of the grid. Figure 1 
illustrates the grid domain and boundary placement in the 
horizontal. Boundaries in the vertical are similarly 


Placed, with bottom and top boundaries at 


, 2243 
aes one 


and 


(KM-1) + (KM-2) 
2 


respectively, where KM is the number of grid points in the 


wemtical. 
i=l i=IM 
j=JM j=IM 
——" <8 North boundary So 
pone --  - - - - e------- ----- 
® o | fs ® 
| | 
| | 
* | e | ° 
| | 
West} boundary East! boundary 
| 
| | 
e ° | e lee © 
® ahs & ies 
pa Rs a 1 
nn South boundary ee 4 we 
i=l i=IM 
oe jel 


Pigugem sms Horizontal grid 


Periodicity in the east-west direction ensures cyclic 


continuity of all parameters. This condition is imposed by 


IS, 








applying the relationships 
ee These ke 


"2,418 9 Sane, ae 


S = 
IM-1,j,k a ee 


and Samy yk = S4,j,k ; 
where 8S = u, v, w, @ and P. Since values of S at points 
in the outer columns are required in solving for the pressure 
term and in calculating divergence, a four point periodicity 
Constraint is necessary. 

Mimorcernce, Pprecvene the £low OF mass and energy through 


the northern and southern boundaries, a condition of no flux 


ls imposed on the v component of the wind. This is accom- 
plished by 
iO ee | Soe 
and 
“i7cM=),k 9 -i,J3M-2,k ° 


Ge Values of v on the outer row are obtained by 
Vi,l,k~ “4,4,k 


and 


ar Vv. 
1,J5M,K 1,JM-3,k . 
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The u and w wind components and @ are projected 


Outward by 

weer a 2,k  «t,3,k 
and See M1, %  “i,gM-2,k * 
oe S = uw, @. 


Values of pressure on the outer rows are obtained by using 


The geostrophic relationship 


op 
Oy 


= —-fu © 


The upper and lower boundaries are treated in the same 
Manner as the north and south boundaries, with the no-£lux 


Condition being imposed on w by 


= Ws ; : 
and i,j, KM thy J) g tls 


Valles Ot swe and Vv on the top and bottom levels are 


obtained by 
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and 


i,j, KM iy Gee 


where S= u,Vv. 


ies nydrostatic relationship, 


g@ 


Zz 0, 


Q/ 
© 


f 


O/ 


Sea 2 ee 


is used to determine values of @ and @® above the top and 


below the bottom boundaries. The application of proper 


boundary conditions was found to be critical. Appendices C 


and D show the development of these conditions in detail. 


Ze 





Ven LNT TTAL CONDITIONS 


Initialization of the u,v,w and @ fields was 
Pecomplished according to the experiment being prepared. 
During the check~out phase of programming the model, a 
Steady-state, geostrophic, zonal flow was used to provide a 
Check on computational stability. Subsequent tests were 
initialized to include gravity waves, both stable and 


unstable. 


A. ZONAL FLOW CASE 

Healey Or the balanced State was tested by inserting 
a lO m sect wind throughout the u field. Initial v 
and w fields were set to zero. Standard atmospheric 
Values of pressure were introduced and, by using an adia- 


batic lapse rate of temperature, the equation 


provided the initial @ field. Since the relaxation 
technique used to obtain the initial ® field proved to be 
extremely efficient, choosing the initial guess ® field as 


iieavoltbhary Constant was convenient. 


Becnay biy WAVE CASES ; 
Linearized gravity wave solutions, accurate for small 
Scale perturbations only, were used to test the model. For 


the stable gravity wave the following relationships describe 


the initial conditions. 
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Am W 


u oy sin Ax cos py cos mz 

Vv ao ee cos Ax sin Ply cos mz 
N +p 

w = Wocos Axcos Py sin mz 

QV = 1S oo sin Ax cos py sin mz. 


Poese speed may be Calculated from the relationship 


Ne + pe? 


ee We ae (15) 


For the gravity wave associated with unstable stratifi- 


Sammon, the initializing relationships are 


_ A NA \ 
ee an a ee) W sin Ax cos fly cos mz 
2 
r N | 
—_ A—F + <a 
V = IL oe W cos Ax sin Ply cos mz 
w = W cos Ax cos pty sin mz 


6 “a — cos Ax cos py sin mz , 
Zz 


where 


@) 


3 
Oz 


| 


S g 
N? = — 
0 
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and 


Dene 2 
2 yee (16) 


e+ pe + m2 


where n is the growth rate of the disturbance. 
mics q@ tield may be initialized for the gravity waves 


Seoetly aS in the zonal flow case. 


ZS 





Vile hbo ULTS 


The balanced state introduced for the one dimensional 
flow case as a check on the computational stability and 
forecast capability of the model should remain indefinitely, 
Subject to the validity of the approximations in the basic 
equations, round-off error in the computer, accuracy of the 
peeunotions at the boundaries, and stability requirements. 
The use of proper boundary conditions was the most important 
factor in determining the Leng Ehnwonetic forecast produced. 
Until the combination of boundary conditions shown in 
waeutonm LY and Appendix C was utilized, the forecasts were 
Varerae tor no longer than 5 minutes. Longest forecasts of 
Che balanced flow were obtained by applying the boundary 
Sendttions that are periodic in x, geostrophic in y, and 
icresctatic in 2Z. 

Accuracy of forecasts was monitored through the calcula- 
Cion of the three dimensional divergence, which should have 
remained very small throughout the prognosis. The smallest 
possible relaxation tolerance consistent with single preci- 
Sion mode capability of the computer was found to allow the 
least amount of divergence. When larger tolerances were 
allowed, a substantial amount of divergence developed on the 
north and south boundaries at mid-levels and increased with 
Hee mies 

A test was conducted to determine the optimum relaxation 


coefficient for the grid interval and grid domain which the 
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model employs. The value selected did not appear too 
Critical, as long as it was between 1.5 and 1.7. 

Poceene Liboe Cxperiment that produced a relatively long 
prognosis, an adiabatic lapse rate of temperature and 
standard atmospheric values of pressure were used to ini- 
tialize the potential temperature field. This produced a 
Slightly unstable average lapse rate of potential tempera- 
ture, with two thin stable layers included near the bottom. 
With these conditions the model forecast for 47 minutes. 

The time step was then decreased from 1 minute to 40 
seconds, with the result that the forecast proceeded for 72 
time steps, which was essentially the same length as the 
previous forecast. In an effort to improve the initial 
balance, pressure values which would lead to a nearly neutral 
lapse rate of potential temperature were inserted in the 
first experiment. The forecast continued for 56 minutes. 
Returning again to the conditions used in experiment one, 
the lapse rate of potential temperature was modified to 
include stable layers across the top and bottom boundaries. 
With these conditions the model produced a forecast for 69 
minutes. 

Introduction of a small perturbation in the form of a 
gravity wave provided another method of observing the model 
behavior. The gravity wave in stable stratification should 
propagate in the x direction only, without growing. With 
unstable stratification the wave should not propagate, and 


should grow at a rate given by equation (16). 
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The wave inserted for the stable case was forecast for 
30 minutes, after which a comparison was made between the 
initial and forecast values of u, v, wand @. With this 
comparison it was possible to observe the amount of spurious 
Growth and the speed of propagation. The relationship of 
Ege initial and forecast wave forms in x for the line j=8, 
k=6 is shown in Figures 2a, 2b, 2c and 2d. The exact solu- 
tion of equation (15) indicated that the phase speed for the 


Case shown was 10.2 m sec, Phase speeds calculated from 


the plotted values were 10.5 m sect, 10.4 msec™! , 


S25 m ee and 10.3 m sect 


for u, v, wand @ respec- 
tively. The difference in phase speeds between w and @ 
Was in no doubt related to the fact that the wave grew. The 
fchernorms in y illustrated in Figures 3a, 3b, 3c and 3d 
show that the wave grew in the y direction but did not 
propagate. Figures 4a, 4b, 4c and 4d show that no propaga- 
mom Of the wave occurred in the 2 direction, but that the 
growth was not smooth. 

hOreecachmCcomolmcadtlon OL initial Conditions applied to 
the wave, growth was apparent. By modifying the lapse rate 
of potential temperature toward adiabatic conditions, the 
growth decreased. 

The gravity wave with unstable stratification was 
compared in the same manner as was the stable case. Figures 
Doe oeC ane 5G show that the wave did not propagate. 
Furthermore, the growth that was experienced by the wave was 


very close to that which was expected. The amplification 
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factor based on the computed growth rate of .1376 x 1074 
sec”’+ was determined to be 11.9 at 30 minutes, and amplifi- 
Saeron £actors of 11.8, 12.0, 12.4 and 11.3 were noted for 


u, v, wand @ respectively. 
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Fig. 5b - The v component of the wave form in the x 
direction for i=S, k=6, unstable case 
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Fig. 5d - The 9 wave form in the x direction for 
j=8, k=G6, unstable case 
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VII. SUMMARY AND CONCLUSIONS 


eh emo Lit meive equation model was applied to 
eemesoscale Grid for two test situations. First, a one 
dimensional, balanced, steady state current was initialized 
to check the computational stability and finite difference 
imneOn the equations Of motion, Then an initial gravity 
Rave perturbation was introduced for both stable and 
emicicaple stratification. 

mnie stable stratification case yielded a reasonably 
pecirdte phase speed Of propagation but some undesirable 
growth of the disturbance was observed, as was discussed in 
mney Previous Section. 

The unSteble stratification case produced nearly correct 
growth rates but some difference in the w and @ growth 
rates may be an indication that a problem exists similar to 
that in the stable case, where w and §@ waves were propa- 
gated at slightly different speeds. Apparently, some small 
error still exists in the program. 

The fact that the growth of the propagating gravity 
wave can be decreased by inserting a more nearly adiabatic 
lapse rate of potential temperature can be explained by 
observing that the buoyancy term in the Boussinesgq equations 
becomes increasingly inaccurate with increased height. By 


assigning to 6. a value of potential temperature from a 


O 
middle level, the error introduced in the buoyancy term is 


minimized. 
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Tice ete boOuUssINescmeduletions are accurate only for a 
Shallow atmosphere, the model could be expected to forecast 
epemaccuLatce byw. the vdeprm were decreased. In such a case 
the problem with the buoyancy term would be further mini- 
mized. 

Further studies using this model should have a provi- 
Seen Lor including vertical compressibility. Additionally, 
Since a primary source of the energy in a thunderstorm—-type 
round US latent heat, application to cumulus scale convec— 
tive activity would be incomplete without the inclusion of 


moisture, 
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APPENDIX A 


EXTRAPOLATED LIEBMANN RELAXATION SCHEME FOR 
THREE DiMbhNoi@ONnaAl, APPLICATION 


The equation to be solved by relaxation is 
a 
Veet = (A-1) 
where 
g 9306 


0 ee 
ys [3 Vs | rn ary 


Writing (A-1) in finite difference form and expanding yields 


2,4 .k ~ ZW Boe ag ee 
(2 Ax)2 


Pi j+2,k Sari, k SOR faa ae 
(2 Ay) 2 


, Padme Paget Pai | 
(2Az)2 Ly J) pes 


Rearrangang terms leads to 
= 2 2 
noe | (2 Ay) 2 Az) (CRD ea ie aL Pi -2,5,x! 
+ (2Ax)* (2A2)* ($4 542.5 + O45, 3-2,4) 


H(2 Ax) -(2 Ay) (Pi 5, K+2 u ee ao) 
- (2 Ax)?(2 Ay) 2(2Az) 2 Pi 3,34 A 
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where A is the fraction 


L 


2 | (24v)7 (2.82)? f (2 Ax) (2Az)- + eae Ay)" 


ond 
Adding and subtracting (Oe aa produces 
sds 
new old 
p = +a | (2Ay)2(2Az2z)°(%.,, .,+9...,. 
iisk  4ydek (2 Ay) “(242)°( Pisa jx + Piz, 5,%) 


+ (2Ax)7(2Az)2( 6, 0 tO, od 


+ (2 Ax)*(2 Ay) ?(, mea 1 Pa, 5, Ku2! 


EN old 
ee OS 2) 2 2? 
A aja 


For successive over relaxation the coefficient ad is 


included, which results in the final form 


new old 


Miedo esa), Piij.xt i | (2Av)2(2.A2)70 i495 


: ele2 jal 


+ (2 Ax) 2(2 Az) *( Pi, 5+2,k rie ony 


+ (2 Ax)7(2Ay) 2( Pi,j,kt2 +Pi,j,x-2) 
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Note that if Ax = Ay = Az, equation (A=-2) reduces to the 


more familiar 


new old 
p a _ da 
i,j,k (1 -@) Doe ae eon at eo 4 
mmm PG ta + Pe. 4 uno 


Zz 
(2 Ax) Fs jan : 
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APPENDIX B 


DEVELOPMENT OF THE FORCING FUNCTION 


iimoOcdCMuromictMttiatnmecenalcuent finite differencing 
Beeween the prognoseic equations and the forcing function of 
the Liebmann relaxation equation, a procedure which was 
recommended by Dr. R. T. Williams (personal communication) 
is used to develop the forcing function. The process is 
Gack Led Out in 6 Steps for each point in the grid. The 
6 steps required to build Ps GK PneCmEeone ing sEunCcIOn 


at i=5, j=5, k=5, are 


i ees oye 


SES) a2 Sas 


4 Gluj6 5.5 3 
a7 2 2 4X SoD 


5 _ 6W)5,6,5 4 
SP els. 2A De 


6 CS GGG a & 
5 5 25 Pez, DAS ‘ 


where 


alt 





g@ 
ehh) = a8 renpots - 
(w) 6, sree, V3 W, 
G(v) = — Lv) - £u + VY V5. ve 
and G(u) = —~ Lu) + Ev + Vv v5. ie 


ae 
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BOUNDARY CONDITIONS AT THE WALLS 


HOcetne sp eome aescrioed by j—-S, the equations o£ motion 


may be written 


du op 
cae ee ers eC u). i eereeee = 
a8 ao ( JA se 5c S%; 3) + NIE a (e 1) 
av = ) oP 
a = WN a nas = = fui 3,k (C~2) 

Theor Ie nly, Sara < 

9 @. 3.k 

ow =- 2p 5, -2f eee (C23) 
Om | Oz 6 

eS kk Ale oan Oo 


Forming the time derivative of divergence from (C-1), (C-2), 
and (C-3) and using finite differencing yields 


au _ 0 u 
ot Se 


nes el MS a TLRS) A Ie 


oo (c-4) 





e conditio® 
(C-5) 






pntroducing ch 
V5 2,k =V¥, 3,5 
which s;mpLLes 
ON — ov 
g 
ine oe 
4, 20% ee 
(c-4) pecomes 
roe _ oe 
: Se ati. 3s% St inde 3s® 
28% 
gon , & 
at at 
3, 40% 4,230 
+ 
2 AY 
Ow . ow 
Se 4 See Ae A Bolhee 
aH (c-6) 
2 
Expansion of (c-6) yields 
yo ao a Tugs O, 4,4. O4,21% 
; WU = WU ' + P 
2 &% Dey 
£( Vqele?: S inet 


e(w) 4,3,%24 
f f ae 


e w) 4,3, tt ” 
2 b@ 
1,2! g 3k? 
: (c-7) 


a U4 ,3, 
* 2 
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Applying the conditions 


ores k (C-8) 


ll 

! 
N 
> 

MK 
Hh 

S 


aay t (C9) 


and 


(Ga10) 


lrcoauces 


2 .. 
V ? ty Ops i. eae ak ‘ 


which is the equation to be solved by relaxation. 


Thus, the conditions for the southern boundary must 
include the relationships given by (C-5), (C-8), (C-9) and 
(C-10). 

A Similar approach produces corresponding conditions 


at the northern boundary. These conditions are 


Leas = eM —o kk 4 (C-11) 
U3 dM-1,k = “i,gM-2,k *¢ (C-12) 
%. IM, OSes eee ae = OS CSAE UM siya (C-13) 
and 
Cn Te = PP; sM-3,k ee owe “i IM-2,k (C-14) 
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For the plane described by k=3, the same process can 
be carried out to develop boundary conditions at the top and 


bottom, which are 








Warr dee _ Skee (C=) 
2Azg 

? 3 ~ eat O96 0; 4,2 t (C-16) 
2Azg 

ae Pe ae jas (C-17) 

Wi,j,KM-1 ~ ““i,j,KM-2 ° (C-18) 

_ 2Azg 6 
i,j,KM i,j, KM-2 8. i,j,KM-1 (e=%2) 
and 
= _ 242g 9 
by Ig Ul i,j, KM-3 0, iL, j,KM-2 = 
(C-20) 
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